package de.unibi.comet.tools;

import java.util.List;

import de.unibi.comet.ac.AbelianPattern;
import de.unibi.comet.ac.AhoCorasickAutomaton;
import de.unibi.comet.ac.StringSetProvider;
import de.unibi.comet.fa.Alphabet;
import de.unibi.comet.fa.CDFA;
import de.unibi.comet.fa.DnaMotifParser;
import de.unibi.comet.fa.MarkovAdditiveChain;
import de.unibi.comet.util.Log;

public class abelian_distribution {

	public static void usage() {
		System.out.println("usage: abelian_distribution <cg-content> <length> <max_mass> <pattern>");
		System.out.println("   <cg-content> cg-content to be used as i.i.d. model");
		System.exit(1);
	}

	public static void main(String[] args) {
		if (args.length!=4) usage();

		double cgContent = Double.parseDouble(args[0]);
		if ((cgContent>=1.0) || (cgContent<=0.0)) {
			System.out.println("invalid cg content");
			System.exit(1);
		}

		
		long steps = Long.parseLong(args[1]);
		long maxMass = Long.parseLong(args[2]);
		
		String pattern = args[3];
		
		Log.getInstance().setTimingActive(true);
		// Log.getInstance().setLogLevel(Log.Level.INSANE);
		Log.getInstance().setLogLevel(Log.Level.EVERYTHING);

		Alphabet alphabet = DnaMotifParser.getDnaAlphabet();
		double[] charDist = new double[4];
		charDist[alphabet.getIndex('a')] = (1.0-cgContent)*0.5;
		charDist[alphabet.getIndex('c')] = cgContent*0.5;
		charDist[alphabet.getIndex('g')] = cgContent*0.5;
		charDist[alphabet.getIndex('t')] = (1.0-cgContent)*0.5;

		// maximal value for product maxMass*states*steps
		long limit = 500000000L;
		
		AbelianPattern ap = new AbelianPattern(pattern);
		long stringCount = ap.countAllStrings();
		
		long nodes = -1;
		long nodesMinimized = -1;
		double timeMinimization = -1.0;
		double timeConstruction = -1.0;
		double timeMacEvolution = -1.0;
		double timeStringGeneration = -1.0;
		double pValue = -1.0;
		
		Log.getInstance().startTimer();
		List<String> l2 = ap.generateAllPatterns();
		Log.getInstance().stopTimer("generate strings");
		timeStringGeneration = Log.getInstance().getLastPeriodCpu();
		
		System.out.println(String.format("abelian pattern: %s --> %d strings", ap.toString(), l2.size()));
		
		StringSetProvider ssp = new StringSetProvider(l2);
		AhoCorasickAutomaton aca = new AhoCorasickAutomaton();
		timeConstruction = Log.getInstance().getLastPeriodCpu();
		aca.build(ssp);
		CDFA cdfa = aca.createCDFA(alphabet);
		nodes = cdfa.getStateCount();
		cdfa = cdfa.minimizeHopcroft();
		timeMinimization = Log.getInstance().getLastPeriodCpu();
		nodesMinimized = cdfa.getStateCount();

		long prod = steps*maxMass*nodesMinimized;
		System.out.println(String.format("product: %d",prod));
		if (prod>limit) {
			System.out.println("steps*states*maxMass too large!");
			System.out.println(String.format("States: %d, maxMass: %d, product: %d", nodesMinimized, maxMass, prod));
			System.exit(1);
		}
		
		MarkovAdditiveChain mac = cdfa.createMAC(charDist);
		double[][] dist = new double[cdfa.getStateCount()][(int)maxMass+1];
		dist[0][0]=1.0;
		mac.setDistribution(dist);	
		
		mac.step((int)steps);
		double[] massDist = mac.getMassDistribution();
		timeMacEvolution = Log.getInstance().getLastPeriodCpu();
		pValue=massDist[(int)maxMass];
		
		StringBuffer sb = new StringBuffer();
		sb.append(">>");
		for (int mass=0; mass<=maxMass; ++mass) {
			sb.append(String.format(" %e", massDist[mass]));
		}
		System.out.println(sb.toString());
		
		System.out.println(String.format("> %s %d %d %d %e %e %e %e %e", pattern, nodes, nodesMinimized, stringCount, pValue, timeMacEvolution, timeStringGeneration, timeConstruction, timeMinimization));
	}
}
